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We investigate the quantum- mechanical localization of X H and 2 H isotopes in the symmetric 
low-barrier hydrogen-bonds of potassium dihydrogen phosphate (KDP) crystals in the paraelectric 
phase. The spatial density distributions of these hydrogen atoms are suspected to be responsible 
for the surprisingly large isotope effect observed for the ferroelectric phase transition in KDP. We 
employ ab initio path integral molecular dynamics simulations to obtain the nuclear real-space 
and momentum-space densities n(R) and n(k) of 1 H and 2 H, which are compared to experimental 
Neutron Compton Scattering data. 

Our results suggest a qualitative difference in the nature of the paraelectic phase in KDP between 
the two isotopes. We are able to discriminate between real quantum derealization and vibration- 
assisted hopping and thus provide evidence for two distinct mechanisms of the ferroelectric phase 
transition in this class of materials. 
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Potassium dihydrogen phosphate (KDP) is the proto- 
type of a wide class of hydrogen-bonded ferroelectrics 
that find use in non-linear optics applications Its 
phosphate units are interconnected by a network of 
hydrogen-bonds (H-bonds). The ferroelectric to para- 
electric phase transition in these ionic crystals is asso- 
ciated with a shift from asymmetric OH- • • O H-bonds 
(fig. [T|) to a symmetric configuration. 

Earlier theories for the phase transition associated the 
symmetry of the paraelectric phase with a classical dis- 
ordering of protons over the two donor sites in each H- 
bond However, the observation of a very large iso- 

tope effect (T C [ 1 H]=122K to T C [ 2 H]=230K 0), strongly 
suggested the presence of quantum effects. As a result, 
the observed disorder was attributed H solely to X H- 
tunneling between the two equivalent donor sites in each 
H-bond. The increased mass due to deuteration reduces 
the tunnel-splitting. As the asymmetric OH- • • O bond 
is favored by electrostatic effects, it would then domi- 
nate the tunneling derealization over a larger range of 
temperatures and lead to a higher T c . 

This tunneling model assumed that the mass effect is 
dominant and that the effective potential energy surface 
(PES) for the 1 H/ 2 H motion remains the same. However, 
pressure-dependent studies on KDP [(j and related sys- 
tems have shown [7] that the tunneling hypothesis is in- 
sufficient to completely account for the observed increase 
in T c . An additional Ubbelohde [H effect is assumed to 
contribute significantly to the increase in T c . 

In this Letter we probe the nature of the paraelec- 
tric phases of protonated and deuterated KDP using 
quantum-mechanical density distributions of the hydro- 
gen nuclei both in real and reciprocal (momentum-) 
space. We use a recently developed ab ini tio p ath- integral 
molecular dynamics (PIMD) method [H, [Toj], which per- 




FIG. 1: Left : The paraelectric KDP unit cell along the polar 
c-axis. The H-bonds between phosphates lie along the a and b- 
axes, respectively. Right : two phosphate molecules with their 
potassium counterions which form the basic building block of 
a KDP unit cell. The ferroelectric polarization is caused by 
a symmetry breaking in the K - P - K distances along the 
c-axis. Bottom panel with a single OH- • • O bond illustrates 
the proton-localization coordinate 5 and the acceptor-donor 
distance 2R (see text). 



mits us to model simultaneously electronic and nuclear 
quantum effects at realistic temperatures and in the ac- 
tual periodic crystal structure. We present quantitative 
results for the localization and derealization of protons 
and deuterons, and thus distinguish between quantum- 
mechanical derealization (tunneling) and classical hop- 
ping of localized particles, as a function of tempera- 
ture. We show that this essential mechanistic differ- 
ence between the ferroelectric transitions in protonated 
and deuterated KDP is key to understanding the huge 
isotope-effect observed in this sytem. 

Solving the nuclear Schrodinger equation on a density- 
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functional theory (DFT) based PES requires prohibitive 
computational resources except for very low-dimensional 
systems, i.e. with few nuclear degrees of freedom. Our 
calculations are based on the path integral technique, 
which represents a convenient way to map the nuclear 
quantum problem to a statistical ensemble average of 
a set of non-quantum replica of the same system. The 
general path-integral method is described in detail else- 
where (nUTs). 

In the common path integral formulation of a quan- 
tum system of n nuclei, the canonical partition function 
Z = Tr exp — f31-L = Tr /3[ T ] is written as an integral 
over discretized paths {R^R^ . . . R/ p ); e 3? 3n } in 
coordinate-space. 
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where /3" 1 = k B T and P E M is a large integer, 
ARCp) = R^ - R(p +1 ) and R( p+1 ) = Rt 1 ). This par- 
tition function is equivalent to that of n (classical) ring 
polymers (one per particle), each composed of P elements 
(beads) connected by classical springs with spring con- 
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The ensemble of n particles for a given 



bead is called a replica. The particles within the p th 
replica interact through the potential V (R^ p ^), which in 
our case is the total electronic energy. 

Every replica is a classical representation to the quan- 
tum system under study, while the ensemble of the P 
coupled replica is isomorphic to the quantum system of 
n atoms. The isomorphism is exploited to compute the 
solution of the n-particle quantum problem via a thermo- 
dynamical phase space sampling of the n coupled classi- 
cal ring polymers. This sampling can be done by means 
of Monte Carlo techniques or molecular dynamics (MD) 
simulations, our choice being the latter (using the CPMD 
code [18]). Expectation values of all quantum observables 
which are local in direct space can be expressed via the 
statistical density of beads, in particular the probability 
density 
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Properties which are represented in the space of the 
momenta p = ftk, however, require the momentum den- 
sity n nuc (k) = Tr (|k)(k| P[t])i an d thus a modifica- 
tion to the original path-integral formulation. Since |k) 
is an eigenstate of the kinetic energy operator T and 
(R|k) =e- ikR , we obtain 

n nuc (k) = e-^ k ) J d 3n R^ d 3n R^ x 

e -ik(RW-R<*>) /( R (1) 5 R(2)). (4) 

with the off-diagonal density matrix element p(R, R'). 
The latter can be sampled within the classical path- 



integral isomorphism via a open polymer, where the 
classical spring between replica R^ 1 ) and R^ 2 ) has been 
eliminated jl6j |. Using this open path integral formal- 
ism at P = 32 within a gradient-corrected DFT frame- 
work pjl [20| we have computed the direct- and recip- 
rocal space density distributions n(r) and n(k) of the 
hydrogen atoms in KDP and DKDP crystals by open- 
ing the paths for only these atoms [24[. The simulation 
temperatures were chosen at least about 20K above the 
respective transition temperatures (140K and 300K), in 
order to compensate for the known temperature issues in 
DFT-based MD calculations [Tcj. Our supercell for the 
paraelectric phases contains four independent phosphate 
units on body-centered tetragonal sites, which gives us 
sufficient statistics (8 protons per cell) but ignores cou- 
pling to phonons away from the zone-center. The lattice 
parameters for each isotope were taken to be the experi- 
mentally measured values in the paraelectric phases close 
to the T c . For KDP a = 7.4264A,c = 6. 931 A, and for 
DKDP a = 7.459A, c = 6.957. Note that the difference in 
volume is only about 1.3%. The total simulated time was 
about lOps (at a O.lfs time-step) ensuring good statistics 
covering both localized and delocalized paths. 

As illustrated in fig. [TJ the KDP crystal contains a 
two-dimensional network of OH- • • O hydrogen bonds be- 
tween the phosphate units, orthogonal to the c-axis along 
which the electric polarization is oriented in the ferroelec- 
tric phase. We characterize the H-bonds by the oxygen- 
oxygen distance 2R and a hydrogen-localization coordi- 
nate S = R — HoiH ' R>OiO u which quantifies the asym- 
metry of the OH- • • O bond. Here, the subscripts l(u) 
stand for the "lower" ("upper") oxygen, Ho t H refers to 
the instantaneous vector from the lower oxygen to the 
hydrogen, and Ro ; o u is the unit vector pointing from 
the lower to the upper oxygen. Each MD step provides a 
set of P values for 5, whose statistical distribution n(S) 
over the PIMD run defines a hydrogen-localization func- 
tion (HLF). The HLF represents the projection of the 
real-space density of the hydrogens along the OH- • • O 
bond. A negative (positive) value of S implies that the 
hydrogen is closer to the lower (upper) oxygen O^ (O u ) 
of the H-bond (see fig. [I]). 

Fig. [2ji shows the HLF for protons in KDP resulting 
from the PIMD simulations. The broad HLF can result 
from either quantum derealization or a classical hopping 
across a barrier in the H-bond. However, a conventional 
Car-Parrinello (CP) simulation at the same conditions 
showed no hopping at all (see the inset in fig.[2k) but, in- 
stead, resulted in a symmetry-broken ferroelectric state. 
Hence, in the time-scale of this simulation, a classical 
proton cannot overcome the OH- • • O barrier at 140 K. 
This implies that the broad HLFs are due to quantum 
derealization of the protons in the H-bonds, leading to 
the symmetric paraelectric phase. 

The computed HLFs of the deuterons in paraelectric 
DKDP at 300 K is shown in fig. Hi (bottom left). In con- 
trast to KDP, we see a diversification of the HLF ranging 
from well-localized off-centered deuterons to considerably 
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FIG. 2: (a) Hydrogen localization functions HLF (left) and 
simulations for 1 H-KDP at T=140K (top) and 2 H-KDP at 
from a conventional MD simulation of the same system; (b) 
(NCS) O. 



2K-5 correlation (right) from 12ps of ab-initio path-integral MD 
300K (bottom) (paraelectric phases). The inset shows the HLF 
Average momentum distributions - computed and experimental 



delocalized ones. An overall average in them would lead 
to a symmetric distribution. Unlike in KDP, this sym- 
metry results from classical disordering of the deuterons 
between the two OH- • • O equilibrium sites. A comple- 
mentary fingerprint of the paraelectric KDP (DKDP), is 
the X H( 2 H) nuclear momentum distribution n(k) along 
the OH- • • O bond from our linear-polymer path-integral 
approach. These are shown in fig. [2]d along with the 
corresponding experimental data from neutron Compton 
scattering experiments [I?}- The good agreement with 
experiment confirms the 1 H-tunneling picture in KDP. 
The n(k) of the deuterons is broader than that of the pro- 
tons, which corresponds to a tighter localization in real- 
space. In fact, the experimental np#](k) corresponds to 
an even higher degree of spatial localization further vali- 
dating the classical disordering picture. 

We have further computed a = yj (& 2 ) _1 as a mea- 
sure of the instantaneous real-space hydrogen localiza- 
tion (see the a scale-bar in fig. Et). While the deuterons 
are more localized than the protons, the spread of either 
species lies between the extremes of a free particle and a 
simple harmonic oscillator. Interestingly, deuteration at 
the temperature and volume conditions of the paraelec- 
tric KDP simulation resulted in only partial localization 
(DKDP at 140 K). The localization is enhanced upon in- 
creasing the temperature. This enhancement is closely 
linked to the so-called geometry-effect of isotopic substi- 
tution and, in fact, a key feature of the mechanism of the 
ferroelectric phase transition in DKDP. 

The geometry effect manifests itself in experiments as 
an increased OH- • • O distance (2R in fig. [I]) in DKDP 
with respect to KDP in the paraelectric phase [6]. Our 
PIMD simulations reproduced this experimental trend 
with an average 2R of 2.49 A in KDP and 2.60 A in 



DKDP. Furthermore, we observed a correlation between 
the instantaneous 2R values and the localization 5 of the 
hydrogens, shown in fig. [2^. At shorter 2R the hydrogen 
prefers the center of the H-bond while for larger 2R the 
most-probable positions of the hydrogen are off-center. 
In fact, for very short H-bonds (2R < 2.45 A) the off- 
center peaks are absent suggesting the disappearance of a 
symmetry-broken phase at high pressures [211 . Thus, the 
probability of finding the hydrogen at the center of the 
H-bond is effectively modulated by the acceptor-donor 
distance 2i?, which determines the tunneling barrier [22|, 
[23[ . The actual tunneling regime corresponds to a rather 
thin range of 2R values where the S coordinate has a 
distinctively larger spread about the center. This range 
(occuring around 2.49 A) is wider for KDP than DKDP 
due to the higher mass of the 2 H. The stabilization via 
tunnel-splitting in KDP is strong enough to soften of the 
ferroelectric mode (that is coupled to the 2R distances). 
Thus, the phase transition in KDP can also be thought 
of as a soft-mode driven transition. This effect is absent 
in DKDP as the higher mass of the 2 H diminishes the 
tunneling probability and hence the tunnel-splitting. 

As the deuterons in the paraelectric phase of DKDP 
are localized, they must be (classically) disordered over 
the two OH- • • O sites. The distinction betwen a frus- 
trated antiferroelectric or dynamically disordered char- 
acter of this state can be made using the individual mo- 
mentum distributions as well as the corresponding real- 
space localization functions (see fig. [3]). The real-space 
picture shows that there are both localized and delocal- 
ized deuterons; the momentum distributions, however, 
all have very similar widths. In some cases, deuterons 
with a broader real-space distribution (e.g. the green 
line in fig. [3^) exhibit an even broader momentum-space 
distribution than more localized ones (e.g. the blue line 
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FIG. 3: (a) Computed momentum and real-space (inset) densities of several deuterons in DKDP at T=300K, and their spreads 
in A (inset); (b) Evolution of the ri]2 H A8) of one specific 2 H in subsequent lps-windows; and (c) the computed spreads (k 2 ) 1 
of n(k) for 1 H (left) and 2 H (right) for a free particle, the harmonic oscillator and KDP at different temperatures. 



in fig.Et). 

This apparent inconsistency is resolved by following 
the the temporal evolution of the real-space density of 
the green hydrogen, shown in fig. [3}}. Initially, the 2 H 
is localized at negative S values (black, red, violet lines), 
but it becomes delocalized after a few picoseconds (blue 
line) and subsequently localizes at positive S (yellow and 
brown lines). At the end of our run, further derealiza- 
tion starts (orange line). This plot illustrates that the 2 H 
is actually localized most of the time, but jumping from 
one oxygen site to the other, passing through a delocal- 
ized quantum state. 

As the 2 H moves towards the bond-center the corre- 
sponding 2R also decreases (see the correlation plot in 
fig. [2^). This not only reduces the barrier to cross the 
bond-center but also leads to a larger real-space width 
in ri[2#](5) [22|. Thus, the HLF is the widest as the 2 H 
passes through the center of the H-bond. Beyond this 
point 2R increases again simultaneously localizing the 
2 H. Hence, the oscillations in the acceptor-donor distance 



2R, effected by rotational/ vibrational modes of the phos- 
phate groups are ultimately responsible for the dynamic 
disordering of the deuterons in paraelectric DKDP. The 
amplitudes of these modes at lower temperatures (in par- 
ticular around the KDP T c of 140 K) are not sufficient to 
enable the hopping process . Thus the ferroelectric phase 
transition in DKDP occurs via a vibration-assisted hop- 
ping of the deuterons across the H-bonds leading to a 
symmetric paraelectric phase. 

In conclusion, we find that the paraelectric phase in 
KDP is stabilized by a derealization state (tunneling) 
of the protons in the OH- • • O hydrogen bonds, whereas 
the paraelectric state of DKDP results from classical dis- 
ordering of the deuterons, leading to a symmetrization 
in a statistical sense. Our results indicate that deuteron 
hopping is strongly coupled to the vibrations of the oxy- 
gens that form the H-bonds (vibration- assisted hopping). 
This fundamental difference in the mechanism for the fer- 
roelectric phase transition turns out to be responsible for 
the huge isotope-effect that is observed experimentally. 
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